Asymptotic behavior of the electron density and the Kohn-Sham potential 
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It is known that the asymptotic decay (\r\ — > oo) of the electron density n(r) outside a molecule is 
informative about its first ionization potential Iq. Contrary to conventional wisdom, we demonstrate 
that this asymptotic decay may not be the same in every direction. When the highest occupied 
Kohn-Sham molecular orbital has a nodal plane, the decay n(|r p | — > oo) on that plane is informative 
about a higher vertical ionization potential U. The implications for density functional theory are 
manifold: the Kohn-Sham potential does not need to have an asymptotic constant shift on the plane 
as previously believed, provided that the second highest occupied orbital energy equals —Ii. The 
effective potential for \J n(r) may diverge asymptotically on the plane. 



In Kohn-Sham (KS) density functional theory (DFT), 
which is perhaps the most widely used technique in elec- 
tronic structure theory, ground-state properties are cal- 
culated via the KS system, consisting of non-interacting 
electrons moving in the local KS potential v s (r) . In prin- 
ciple, the KS potential v s (r) ensures that the electron 
density n(r) of the non-interacting KS system is the same 
as that of the physical, interacting system. Exact prop- 
erties of n{r) and v s (r) have played - and continue to 
play - a crucial role in constructing and improving ap- 
proximations. 

Both the (square root of the) density and the KS or- 
bitals ifjk{r) obey Schrodinger-type equations, 

(-^ v2 +«oxt(r) +v eS (r)\ \/n(r) = -lov^ (1) 

(-^V 2 + u ext (r) + v Hxc (r)\ Mr) = £ki>k(r), (2) 

where v s = v cxt + wh X c and Iq = — Eq is the first 

ionization potential, and J^k IV'/c( r )| 2 = n(r). Here the 
external potential v ex t(f) is a nuclear potential going to 
zero at large distance like —Z/r where Z is the total 
charge of all nuclei and r is the distance from the center 
of nuclear charge. In this case, according to Eqs. ([T])-([2]), 
the asymptotic [r — ¥ oo) decay of ^n(r) and ^fc(r) is 

^ e - v /2(/o+«cff(oo))r and _ e -V 2 (-^+^Hxc(oo))r respec _ 

tively. Both the effective potential v c ff(r) for ^n(r) 
and v s (i.e. the Hartree-exchange-correlation potential 
^Hxc( r ')) have been thought to go to asymptotically. 
In the KS case there is, with the local KS potential, 
no coupling of the eigenvalue equations, as there is 
in the case of Hartree-Fock (HF) pQ. Therefore the 
density decays as the square of the highest-occupied 
molecular orbital (HOMO), leading to the identification 
I = -eh 0| In case the KS HOMO has a nodal 
plane (HNP), a very straightforward argument was 
given by Wu et al. [3] that the KS potential should go 
to a negative constant for asymptotic points r p — > oo in 
that plane. The KS density in that plane is governed 
by the HOMO-1, assuming HOMO-1 docs not have 
the same nodal plane. Wu et al. [3] made the common 



assumption that the exact interacting density has 
the same asymptotic behavior everywhere, and they 
deduced that the decay of the HOMO-1 in the HNP 
mus t be equal to the decay of t he total density, implying 
-yj2(~e H _i +i>Hxc(r p ->• oo))r p = -^/2l^r p , so that 
VHxc( r p oo) should tend to the negative constant 
Iq + eii-i = —(eh — £h—i)- On the other hand, it 
had earlier been argued, and numerical evidence had 
been provided, that the optimized effective potential 
method for the exact exchange model of Kohn-Sham 
(xOEP) leads to an asymptotic constant in the HNP, 
but positive [SHE]. Since it was pointed out that this 
behavior of the Kohn-Sham potential has a significant 
effect on the orbital energies of particularly the higher 
lying unoccupied orbitals, with large consequences for 
the excitation energies calculated with time-dependent 
DFT the matter has practical importance. 

In this Letter we show that: (i) the asymptotic decay 
of the exact interacting density may be different in differ- 
ent directions: if there is a KS HOMO nodal plane, the 
decay «(|t p | — > oo) in that plane will be different (faster) 
than the decay outside the plane. We support our proof 
with results from a simple interacting model that fully 
captures the essence of the problem and is analytically 
solvable; (ii) the potential v e g(r) of Eq. ([!]), while nor- 
mally going to zero asymptotically, exhibits a different 
asymptotic behavior in directions where the density de- 
cays differently: in the KS HNP it can tend asymptot- 
ically to a constant or even to infinity. We then give 
arguments that: (iii) the KS potential can still go every- 
where asymptotically to zero, provided that, if there is a 
HNP, the KS orbital energy of the HOMO-1 is equal to 
the second vertical ionization potential, eu—i = —Ii- 

Asymptotic behavior of the exact density - Although it 
is commonly assumed that the exact density decays ev- 
erywhere in the same way, it should be recognized that 
this is not always true. First of all, it can be easily shown 
for a noninteracting electron system (and also for the 
Hartree-Fock approximation - see below), that if there 
is a nodal plane in the HOMO, the density decays dif- 
ferently in that plane. For the exact density of an in- 
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teracting system, where a configuration mixing will in- 
volve many configurations, the matter is more subtle. To 
study the density decay in the general interacting case, 
we express the exact iV-electron wavefunction and 
the exact density in terms of the Dyson orbitals dt (x) , 



di{x) 
n(x) 



N- 1 ' 2 Y J d t {x)^- 1 (2- ■ ■ N), 

i 

J *f _1 (2 • • • N)*V$(x, 2 • • • N)d2 ■■■dN, 

OO 

El*(*)| a , ( 3 ) 



i=0 



where the 'I'^ -1 are the exact (ion) N — 1 states and 
x = r,s. The sum over i goes over both the spin-f 
and spin-j. Dyson orbitals. If e.g. x — (r,s =t) then 
only the spin-f Dyson orbitals are nonzero at x and con- 
tribute to n(x) = n(r,t) (= (l/2)n(r) in closed shell 
systems). Each state of the ion is associated with a 
one-particle wavefunction, its Dyson orbital. These or- 
bitals constitute a nonorthogonal nonnormal, in general 
linearly dependent set. We define the conditional ampli- 
tude $(2 • • • N; x) and associated quantities, 



$(2---N;x) = 
n cond {x 2 \x) = 
v cond {x) = 



y/n{x)/N ' 
(N-1) J |^(2 ■ - • A^|a;)| 2 cl3 - ■ • cW, 
cond (x 2 \x) 



r2 



-dx 2 



(4) 



$(2 • • • N; x) is a normalized N—l electron wavefunction 
depending parametrically on the position x. Its square 
describes the probability distribution of electrons at po- 
sitions 2 ■ ■ ■ N when one electron is known to be at cc. Its 
associated one-electron density n cond (x2\x) is the density 
of the other electrons at position x 2 when one electron is 
at x, which is the normal one-electron density n(x 2 ) plus 
the full exchange-correlation hole surrounding position x. 
Projecting the Schrodinger equation H'^'&q = Eq^>q 
against , I'^ V ~ 1 (2 • • • N) and using the expansion of Eq. |3| 
one obtains the usual equations for the Dyson orbitals, 



1. 



-V 2 



di(x)+j2* 

k=a 



k (x)d k (x) = -Iidi(x). 



(5) 

Davidson (KD) [9] pointed out 
the coupling integrals Xik(x) = 



Katriel and 
that, due to 

(*r _1 |Ei>il» , i-» , r 1 |*fc r " 1 )2..iv» the exponential 
decay of the coupled Dyson orbitals will be the same. 
This was demonstrated by Handy et al. PQ for the 
analogous case of coupling of the Hartree-Fock orbitals 
by the exchange term. The first Dyson orbital do will 
have exponential decay ~ e _v/27 ° r multiplied by a factor 
r@ with /3 — Z — N + 1, due to the —Z/r decay of w oxt 
and the (N — l)/r decay of the coupling term. KD find 



that higher Dyson orbitals which have nonzero Xio with 
the first Dyson orbital will have decay r^~ L e -^i« r 
with L* > 2. Dyson orbitals that are not connected to 
do will have different exponential decay, governed by 
the eigenvalue of the first orbital in such a connected 
set (which is disjunct from other sets). Considering the 
expansion of the density in Dyson orbitals in Eq. (J3j) , 
KD have noted that, if the density decays for \r\ — > oo as 
the most slowly decaying term |do(a;)| 2 , its exponential 



decay would be ~ e" 



-2V27or 



. Levy, Perdew and Sahni [3] 



proved this exponential decay in a different way, thereby 
proving that the leading term and not the infinite 
summation of Eq. ([3| determines the density decay. The 



\do(r)\< 



-2V2To r 



is considered 



result n(\r\ — ► oo) 
well established. 

This picture changes if the KS HOMO has a nodal 
plane. The common thinking is that still the interact- 
ing density decays everywhere in the same way due to 
correlation effects [4]. Instead, by analyzing the Dyson 
orbitals we can see that a KS HNP also signifies spe- 
cial behavior of the interacting density. A nodal plane 
is typically due to antisymmetry with respect to a sym- 
metry plane (consider, e.g., ethylene or benzene 0E]), 
but a HOMO nodal surface can occur also in more gen- 
eral situations. In the case of a symmetry plane, the 
exact interacting states of the molecule are either sym- 
metric or antisymmetric with respect to the plane. For 
example, the ground state wavefunction corresponding 
to a closed shell configuration is totally symmetric with 
respect to that plane, while the first ion state will be 
antisymmetric (the KS first ion state surely will be so, 
and we will only consider the case that the same holds 
for the exact ion state). For points r p in the HNP the 
conditional amplitude $(2 • • • N\x p ) will be symmetric 
with respect to the plane. Therefore, the matrix ele- 
ment (* v_1 |<i>(2 • • • N\x p )) 2 ,, N will vanish, so that the 
first Dyson orbital is zero in the plane: 

d (x p ) = v /n0r p )(* o v - 1 |$(2 • • • N\x p )) 2 .. N = (6) 

When we consider the asymptotic behavior of higher 
Dyson orbitals, it is clear that with do(x p ) — 0, the cou- 
pling to do in Eq. |5]) for points in the HNP will be zero 
for any higher Dyson orbital di>o- Therefore the decay in 
the HNP of those higher Dyson orbitals is not governed 
by d - We then turn to d\(x p ). If the corresponding 
excited ion state \E'^ V_1 has the same symmetry with re- 
spect to the HNP as , this orbital will not be zero 
in the plane. It will decay in the plane like e -\/27i>p ; 
since in the eigenvalue equation ^ all the terms except 
those from V 2 can be neglected in the asymptotic region 
compared to I\. All higher Dyson orbitals that are con- 
nected to d% in the HNP will decay similarly. Accord- 
ing to Eq. ([3]) the density will, for points in the HNP, 



decay as \di(x p ) 



e 2 v27i r P 5 barring again, as for 



other directions, fortuitous cancellation of this behavior 
by the infinite summation. As an addendum to the re- 
sults of Ref. [I , it may be shown in the same way that in 
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the HNP of a HF atom the other orbitals decay like the 
HOMO— 1. We conclude that the exact density will have 
a different decay in the HNP than in general directions. 
We can prove this if the HNP is a symmetry plane, but 
we expect similar behavior to occur generally in case of 
an asymptotic nodal surface of the HOMO. 

As this result goes against the current belief, we sup- 
port it with a simple fully interacting case that can 
be solved analytically and captures the mathematical 
essence of the problem. Consider two spinless inter- 
acting electrons (with standard Coulomb 1/r interac- 
tion) in an external parabolic potential v ex t(r) = \u 2 r 2 . 
As well known, the corresponding hamiltonian is sepa- 
rable into center-of-mass R = \{r\ + ra) and relative 
r i2 = r 2 ~ r i coordinates, so that its exact wavefunc- 
tion reads (^1,^2) = £(R)</>(t"i2). With spinless elec- 
trons, the spatial wavefunction must satisfy r 2 ) = 
— ^0 (r2,ri), which implies that the ground state corre- 
sponds to the £12 — 1 spherical harmonic for the rel- 
ative vector r\2- We have then 3 degenerate ground- 
state wavef unctions, and we choose one of them by fixing 
m 12 = 0: this way, we obtain an interacting density with 
a symmetry plane like the one encountered in molecules. 
With Uext(f) — ^w 2 r 2 , the way the ionization-energy in- 
formation is embodied into the asymptotic of the den- 
sity is obviously different than the molecular case: we 



with E$ 



E, 



N-l 



have that n(\r\ — » 00) 

w(| + q), where q E K. The values of u for which q 
is integer correspond to analytical solutions of the inter- 
acting N — 2 hamiltonian [TU]. For £12 = 1, ui = \ is 
one of those. We find that the corresponding interacting 
density decays as n{\r\ — > 00) ~ r i e~ ujr everywhere ex- 
cept on the z = plane, where it decays as ~ r 2 e _ " r , 
thus carrying a different ionization-energy information. 
This can be further confirmed by looking at the corre- 
sponding effective potential v c g(r). In fact, as it was 
pointed out by LPS |3], Eq. (TJ) is valid for any binding 
^ext(^)j with v e s(r) expected to go to zero asymptoti- 
cally everywhere (it is only the Coulombic nature of the 
electron-electron repulsion that matters for v c g). Using 
our analytical density we have calculated by inversion 



Ves(r) 



v e xt{ r ) — Iq- The results are shown 



in Fig. [T] we clearly see that v e g(r) goes to zero asymp- 
totically everywhere except on the plane z = 0. In the 
following we analyze in general the behavior of v e g(r) 
close to the HNP. 

The effective potential for yjn - The potential v e g(r) 
of Eq. ([IJ can be written in the form [3l [TT] 



Ve«(x) 



v cond( x) + 1 (Va$(2 . . . jvl^lv^^ • • • N\x)) 



tN-1 



K Ar " 1 l$(2- 



($(2---N\x)\H 

v cond (x)+v km (x) + v N - 1 (x) 



■N\x)) 



(7) 



LPS [3j stressed that each term in Eq. ^ is everywhere 
nonnegative and should tend to zero asymptotically. In 
fact, v cond (x) [see Eq. Q], being the repulsive Coulomb 




FIG. 1: The effective potential v c g(r) for sjn{r) in the case 
of the exactly solvable model described in the text, for larger 
and larger rasa function of 6 — arccos(^). 



potential of a localized charge distribution of (N — 1) 
electrons, decays like (N — l)/r. The third term of v e g, 
jj^ 1 , is positive since in general $ will not be the ground 
state wavefunction of the ion, so its expectation value 
will be larger than Eq' 1 . When r 00 it has been 
inferred that the conditional amplitude collapses to the 
ion ground state [9] (when s =t then $ will collapse 
to the Ms = —1/2 state of the doublet ion), so that 



v N-i., r 



(|r| — > 00) — > 0. The second term, v ' , is mani- 
festly positive and is expected to go to zero asymptoti- 
cally since the derivative of $ with respect to r when the 
reference electron is very far becomes zero ($ remains 
constant - the ion ground state - under small change of 
r at 00). 

LPS [3] and KD [9] carefully note that the above results 
do not hold if ^q 1 is not "accessible" because do (x) = 
for reasons of spin symmetry. However, for points r p in 
the HNP do(x p ) — because of spatial symmetry. By ex- 
panding the conditional amplitude $(2 • • • N\x) in terms 
of the exact N — 1 states, 



t>(2---N\x) = J2 



V n ( x ) 



1 (2-.-iV), (8) 



we see that, since on the HNP do = and n(\r p \ — > 00) 
|rfi(fp)| 2 , on the plane the conditional amplitude tends 



asymptotically to the first-excited ion state, $ 
implying that 



N-l 



.N-l, 



-> 00) = E?- 1 - ES 



= h-h 



(9) 



This is a positive constant appearing in the asymptotics 
of v c s only on the HNP. But this is not the end of the 
story: the second term in Eq. v km , can also be non 
zero at infinity: when crossing the HNP, the asymptotic 



conditional amplitude changes from to ^\ , so 

that the r-derivative of $ perpendicular to the plane can 
be nonzero on the HNP also when \r\ — > 00. Its actual 
value depends on how do{r — > r p ) goes to zero when 
approaching the nodal plane. For example, if we make 
the simple assumption that asymptotically, in spherical 
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coordinates, do ~ /(cos0)e 
di ~ e _v/ ^ r , by writing v km 



m ° r with /(0) = 0, 
in the form [TTHIl] 



and 



|W a (:r)| 2 |Vn(r) 
n(r) 8n(r) 2 



it is easy to see that 



\r p -> oo) 



/27^")r 



(10) 



(11) 



showing that u fcm can go asymptotically to infinity on 
the HNP if /'(0) ^ (as, e.g., in the case of n orbitals). 
A simple illustration of this fact can be obtained by 
calculating analytically the potential v c r{t) corre- 
sponding to the non-interacting density of 3 electrons 
in the configuration ls 2 2p\ in the external potential 
foxt = —\- There, we can clearly see that i> e ff goes 
to zero asymptotically everywhere except on the plane 
z = 0, where it goes to infinity, as predicted by Eq. 
This asymptotic diverging behavior of v e s on the plane 
is perfectly compatible with an analytical, well-behaved 
density. The key point is that when we project Eq. ([I]) on 
the plane, we have to take into account also S7 2 y/n in the 
direction perpendicular to the plane. Usually, this term 
is zero when r — > oo, but when there is a HNP this is not 
the case. In the exact interacting case of Fig. [T] we have, 
on the plane, v km (oo) — const., and we have verified 



analytically that v e s(r p — > oo) = I± — Iq 



.kin 



(oo). 



The Kohn-Sham potential - What happens to the KS 
potential for \r p \ ~ > oo? We have seen that the argument 
of Ref. 4 for a negative constant has disappeared since 
the exact density has a different decay in the HNP than 
elsewhere. Since the KS density in the HNP is domi- 
nated by the HOMO-1, the requirement that the KS sys- 
tem gives the exact density implies that, on the plane, 
tpM-ii^p —> oo) — > e~^ 2 ^ L " r . If v s goes to zero asymp- 
totically everywhere, also on the plane, this would imply 
(-H-i = —Ii- The results of Refs. [5j [6] for a positive 



constant are for the exact-exchange model, and do not 
imply that it also exists for the full KS potential [TH [15] . 
Conversely, it is hard to prove that a positive asymp- 
totic constant C does not exist for v s , since it would be 
compatible with the e -v / 27i> decay of the HOMO— 1 if 
e H -i = -I\ + C (with < C < h — h)- It is inter- 
esting to observe that the few accurate (but not exact) 
calculations that are available for the KS orbital energies 
[131 US] suggest e/f_i = — Jj and thus C — 0, but the 
verification is not very accurate. The KS energies are 
only obtained to ~ 0.05 eV accuracy, and the agreement 
of the first few upper valence orbital energies with the 
first exact primary ionization energies is of the order of 
~ 0.1 eV anyway [13]. The calculations of the "exact" 
orbital energies use KS potentials that reproduce accu- 
rate CI densities. There is no indication that the KS 
potential, which uniquely corresponds to the CI density, 
exhibits any tendency towards a constant C at infinity in 
the HNP. However, such features in the KS potential are 
extremely difficult to obtain with good numerical preci- 
sion. In any case, the order of magnitude of C obtained 
in the exact-exchange model [5J [5] cannot apply to v s in 
view of the accuracy within which eu—i = —I\ for the 
avalaible data. 

Conclusions - We have proved that the exact interact- 
ing density does not necessarily have a uniform asymp- 
totic decay; it carries different ionization-energy informa- 
tion in directions where the KS HOMO has a nodal plane. 
We have investigated the implications for DFT. The ef- 
fective potential for *Jn will in general go to a constant 
in the HNP, or even diverge asymptotically. Our results 
strongly suggest that the exact KS ejj-i = —I\ and that 
the exact v s uniformly tends to zero asymptotically. 
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